Astron. Nachr. / AN 000, No. 00, 1 -0(0000) / DOI please set DOI! 



Equatorial magnetic helicity flux in simulations with different gauges 
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We use direct numerical simulations of forced MHD turbulence with a forcing function that produces two different signs 
of kinetic helicity in the upper and lower parts of the domain. We show that the mean flux of magnetic helicity from 
the small-scale field between the two parts of the domain can be described by a Fickian diffusion law with a diffusion 
coefficient that is approximately independent of the magnetic Reynolds number and about one third of the estimated 
turbulent magnetic diffusivity. The data suggest that the turbulent diffusive magnetic helicity flux can only be expected to 
alleviate catastrophic quenching at Reynolds numbers of more than several thousands. We further calculate the magnetic 
helicity density and its flux in the domain for three different gauges. We consider the Weyl gauge, in which the electrostatic 
potential vanishes, the pseudo-Lorenz gauge, where the speed of light is replaced by the sound speed, and the 'resistive 
gauge' in which the Laplacian of the magnetic vector potential acts as resistive term. We find that, in the statistically steady 
state, the time-averaged magnetic helicity density and the magnetic helicity flux are the same in all three gauges. 
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, ^, ■ 1 Introduction 



The generation of magnetic fields on scales larger than the 
eddy scale of the underlying turbulence in astrophysical 
bodies has posed a major problem. Magnetic helicity is be- 
lieved to play an important role in this process (Branden- 
burg & Subramanian 2005a). The magnetic helicity density, 
defined by A ■ B, where B = V x A is the magnetic 
field and A is the corresponding magnetic vector poten- 
tial, is important because at large scales it is produced in 
many dynamos. This has been demonstrated for dynamos 
based on the a effect (Shukurov et al. 2006, Brandenburg 
et al. 2009), the shear-current effect (Brandenburg & Sub- 
ramanian 2005b), and the incoherent a-shear effect (Bran- 
denburg et al. 2008). The volume integral of the magnetic 
helicity density over periodic domains (as well as domains 
with perfect-conductor boundary conditions or infinite do- 
mains where the magnetic field and the vector potential de- 
cays fast enough at infinity) is a conserved quantity in ideal 
MHD. This conservation is also believed to be recovered 
in the limit of infinite magnetic Reynolds number in non- 
ideal MHD (Berger 1984). This implies that for finite (but 
large) magnetic Reynolds numbers magnetic helicity can 
decay only through microscopic resistivity. This would in 
turn control the saturation time and cycle periods of large- 
scale helical magnetic field which would be too slow to 
explain the observed variations of magnetic fields in astro- 
physical settings, such as for example the 1 1 year variation 
of the large-scale fields during the solar cycle. 
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A possible way out of this deadlock is provided by 
fluxes of magnetic helicity out of the domain (Blackman & 
Field 2000, Kleeorin et al. 2000). In the case of solar dy- 
namo, such a flux could be out of the domain, mediated 
by coronal mass ejections, or it could be across the equa- 
tor, mediated by internal fluxes within the domain. Several 
possible candidates for magnetic helicity fluxes have been 
proposed (Kleeorin & Rogachevskii 1999, Vishniac & Cho 
2001, Subramanian & Brandenburg 2004). 

In this paper we measure the diffusive flux across the 
domain with two different signs of magnetic helicity. This 
measurement however poses an additional difficulty, due 
to the fact that neither the flux nor the magnetic helic- 
ity density remain invariant under the gauge transforma- 
tion A — > A + VA, up to which the vector potential is 
defined. This constitutes a gauge problem. This problem, 
however, does not arise in homogeneous (or nearly homoge- 
neous) domains with periodic or perfect-conductor bound- 
ary conditions, or in infinitely large domains where both the 
magnetic field and vector potential decay fast enough at in- 
finity. In these cases the volume integral of magnetic he- 
licity is gauge-invariant, because surface terms vanish and 
V • B = 0, so that / B ■ V A dV = - J A V • B dV = 0. 
However, in practice we are often interested in finite or open 
domains with more realistic boundary conditions. Also, if 
we are to talk meaningfully about the exchange of magnetic 
helicity between two parts of the domain we need to evalu- 
ate changes in magnetic helicity densities locally even if the 
integral of the magnetic helicity density over the whole do- 
main is gauge-invariant. An important question then is how 
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to calculate this quantity across arbitrary surfaces in numer- 
ical simulations. Ideally one would like to have a gauge- 
invariant description of magnetic helicity. A number of sug- 
gestions have been put forward in the literature (Berger & 
Field 1984,Subramanian& Brandenburg 2006). In practice, 
however, calculating the gauge-invariant volume integral of 
magnetic helicity poses an awkward complication and may 
not be the quantity relevant for dynamo quenching (Sub- 
ramanian & Brandenburg 2006). In this paper, to partially 
address this question, we take an alternative view and try 
to compare and contrast the magnetic helicity and its flux 
across the domain in three different gauges that are often 
used in numerical simulations. 

2 Model and Background 

The setup in this paper is inspired by the recent work of 
Mitra et al. (2009), who considered a wedge-shaped do- 
main encompassing parts of both the southern and north- 
ern hemispheres. Direct numerical simulations (DNS) of 
the compressible MHD equations with an external force 
which injected negative (positive) helicity in the northern 
(southern) hemispheres shows a dynamo with polarity re- 
versals, oscillations and equatorward migration of magnetic 
activity. It was further shown, using mean-field models, 
that such a dynamo is well described by an a 2 dynamo, 
where a has positive (negative) sign in the northern (south- 
ern) hemisphere. However, the mean-field dynamo showed 
catastrophic quenching, i.e., the ratio of magnetic energy to 
the equipartition magnetic energy decreases as i?" 1 , where 
R m is the magnetic Reynolds number. Such catastrophic 
quenching could potentially be alleviated by a mean flux 
of small-scale magnetic helicity across the equator (Bran- 
denburg et al. 2009). Diffusive flux of this kind has pre- 
viously been employed in mean-field models on empirical 
grounds (Covas et al. 1998, Kleeorin et al. 2000). Using a 
one-dimensional mean-field model of an a 2 dynamo with 
positive a in the north and negative in the south, it was 
possible to show that for large enough values of R m catas- 
trophic quenching is indeed alleviated (Brandenburg et al. 
2009). However, three questions still remained: 

1. Can such a diffusive flux result from DNS? 

2. Is it strong enough to alleviate catastrophic quenching? 

3. When is it independent of the gauge chosen? 

In this paper we provide partial answers to these questions. 

We proceed by simplifying our problem further, both 
conceptually and numerically, by considering simulations 
performed in a rectangular Cartesian box with dimensions 
L x x L y x L z . The box is divided into two equal cubes along 
the z direction, with sides L x = L y = L z /2. We shall refer 
to the xy plane at z = as the 'equator' and the regions with 
positive (negative) z as 'north' and 'south' respectively. We 
shall choose the helicity of the external force such that it has 
negative (positive) helicity in the northern (southern) parts 
of the domain. All the sides of the simulation domain are 



chosen to have periodic boundary conditions. The slowest 
resistive decay rate of the mean magnetic field is r\k\, where 
r\ is the microscopic magnetic diffusivity and k\ = n/L z is 
the lowest wavenumber of the domain. 

We employ two different random forcing functions: one 
where the helicity of the forcing function varies sinusoidally 
with z (Model A) and one where it varies linearly with 
z (Model B). This also leads to a corresponding variation 
of the kinetic and small-scale current helicities in the do- 
main. Model A minimizes the possibility of boundary ef- 
fects, while Model B employs the same profile as that used 
in an earlier mean-field model (Brandenburg et al. 2009). 
The typical wavenumber of the forcing function is chosen 
to be fc f = 20fci in Model A and k f = 16fci in Model B. An 
important control parameter of our simulations is the mag- 
netic Reynolds number, R m — u lms /rjk{, which is varied 
between 2 and 68, although we also present a result with a 
larger value of R m . This last simulation may not have run 
long enough and will therefore not be analyzed in detailed. 

We perform DNS of the equations of compressible 
MHD for an isothermal gas with constant sound speed c s , 

D t U = -c 2 Vlnp+^JxB + F visc + f, (1) 

Amp = -V-U, (2) 
d t A = U x B - rinoJ - Vip, (3) 

where F visc = (p/p)(V 2 U + |VV • U) is the viscous 
force when the dynamic viscosity p is constant (Model A), 
and F visc = v(S/ 2 U + | VV • U + 2S In p) is the viscous 
force when the kinematic viscosity v is constant (Model B), 
U is the velocity, J = V x B/ po is the current density, po 
is the vacuum permeability (in the following we measure 
the magnetic field in Alfven units by setting //o = 1 every- 
where), p is the density, ip is the electrostatic potential, and 
D t = dt + U ■ V is the advective derivative. Here f(x,t) is 
an external random white-in-time helical function of space 
and time. The simulations were performed with the PENCIL 
CODlf^, which uses sixth-order explicit finite differences in 
space and third order accurate time stepping method. We 
use a numerical resolution of 128 x 128 x 256 meshpoints. 

These simulations in a Cartesian box capture the essen- 
tial aspects of the simulations of Mitra et al. (2009) in spher- 
ical wedge-shaped domains. In particular, in this case we 
also observe the generation of large-scale magnetic fields 
which show oscillations on dynamical time scales, reversals 
of polarity and equatorward migration, as can be seen from 
the sequence of snapshots in Fig.Q]for a run with R m = 68. 
Here we express time in units of the expected turbulent dif- 
fusion time, T = (rjtokf)^ 1 , where rj t Q — u rms /3fcf is used 
as reference value (Sur et al. 2008). 

Below we shall employ this setup to study the magnetic 
helicity and its flux. We shall discuss the issue of gauge- 
dependence in Sect. [5] 



http://www.nordita.org/software/pencil-code/ 
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Fig. 1 Visualization of the B y component of the magnetic field on the periphery of the domain at different times showing 
the migration of magnetic patterns from the top and bottom boundaries toward the equator. Time is measured in turbulent 
diffusion times, T = (r/tokl) -1 , where rfto = w rms /3fcf is used as reference 



3 Magnetic helicity fluxes 

Let us first summarize the role played by magnetic helic- 
ity and its fluxes in large-scale helical dynamos. The sim- 
plest case is that of a closed domain, i.e., one with periodic 
or perfect-conductor boundary conditions. In the spirit of 
mean-field theory, we define large-scale (or mean) quanti- 
ties, denoted by an overbar, as a horizontal average taken 
over the x and y directions. In addition, we denote a vol- 
ume average by angular brackets, (•}. The magnetic helicity 
density is denoted by 

h M = A ■ B, (4) 
and its integral over a volume V is denoted by 
1 



H M EE (h M ) EE 



V 



A-BdV. 



(5) 



In general the evolution equation of h M can be written down 
using the MHD equations, which yields 

d t h M = -2E ■ B -V • ;F H , (6) 

where 

:F H = E x A + ipB (7) 

is the magnetic helicity flux and E is the electric field, 
which is given by 

E = -UxB + r)J. (8) 

Given that our system is statistically homogeneous in the 
horizontal directions, we consider the evolution equation for 
the horizontally averaged magnetic helicity density, 

-M 



d t h =-2r)J-B-V T , (9) 

where the contribution from the full electromotive force, 
UxB, has dropped out after taking the dot product with B. 
However, the mean electromotive force from the fluctuating 
fields, 6 = u x b, enters the evolution of the mean fields, so 
this contribution does not vanish if we consider separately 

— M 

the contributions to h that result from mean and fluctuat- 
ing fields, i.e. 

d t h™ = 2£-B-2r)J-B- 



dthf = 
where 



j-f = e x a 



25 ■ B - 2f]j ■ b - V ■ 



E x A + ^B, 

iJb, 



(11) 

(12) 
(13) 



and * = * 

In mean-field dynamo theory one solves the evolution 
equation for B, so ;F m is known explicitly from the ac- 

— M 

tual mean fields. However, the evolution equation for h { 
is not automatically obeyed in the usual mean-field treat- 
ment. This is the reason why in the dynamical quenching 
formalism this equation is added as an additional constraint 



T-M 



O-T-M 



equation. The terms /iJT and j b ss k^hf are coupled to 
the mean-field equations through an additional contribution 
to the a effect with a term proportional to kfh { . However, 

the coupling of the flux term is less clear, because there 
are several possibilities and their relative importance is not 
well established. 

In this paper we are primarily interested in T { across 
the equator. We assume that this flux can be written in terms 
of the gradient of the magnetic helicity density via a Fickian 
diffusion law, i.e., 



T i = -K{ V/l f , 



(14) 



(10) 



where K{ is an effective diffusion coefficient for the mag- 
netic helicity density. 

There are several points to note regarding Eq. (fl4l . 
Firstly, both the magnetic helicity and its flux are gauge- 
dependent. Hence this expression should in principle de- 
pend on the gauge we choose. However, as catastrophic 
quenching is a physically observable phenomenon, it should 
not depend on the particular gauge chosen. Secondly, we re- 
call that Eq. (fT~4T > is purely a conjecture at this stage, and it is 
the aim of this paper to test this conjecture. Thirdly, Eq. ( fT~4T > 
is not the only form of flux of magnetic helicity possible. 
Two other obvious candidates are the advective flux and 
the Vishniac-Cho flux (Vishniac & Cho 2001). However, 
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Table 1 



the 



Dependence of B , normalized by B*, 
slopes of the three terms on the RHS of Eq. (fTTb . normalized 
by ?7toi3g q , as well as the value of Kf/»jto- 



Run 




B 2 


2£ B 


2 V j ■ b 


_ -=H 




Bl 


2 


1.1 


9.42 


-9.38 


-0.04 


0.41 


B2 


5 


2.2 


11.18 


-11.14 


-0.04 


0.34 


B3 


15 


2.0 


4.54 


-4.52 


-0.02 


0.27 


B4 


33 


1.7 


2.28 


-2.27 


-0.01 


0.31 


B5 


68 


0.8 


1.15 


-1.12 


-0.03 


0.34 




-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 



none of them can be of importance to the problem at hand, 
because we have neither a large-scale velocity (thus ruling 
out advective flux) nor a large-scale shear (thus ruling out 
Vishniac-Cho flux). 

4 Diffusive flux and R m dependence 

Let us postpone the discussion of the complications arising 
from the choice of gauge until Sect. [5] and use the resistive 
gauge for the results reported in this section, i.e. we set 



ijj = i]V ■ A. 



(15) 



We then calculate J^f and h { as functions of z from our 
simulations, time-average both of them and use Eq. (fT4l to 

calculate Kf from a least-square fit of T i versus — V/i f 
within the range —1.3 < k\Z < 1.3. The values of Kf as a 
function of R m is given in the last column of Table [JJ 

In order to determine the relative importance of equato- 
rial magnetic helicity fluxes, we now consider individually 
the three terms on the RHS of Eq. dTTb . Within the range 
— 1.3 < k\z < 1.3, all three terms vary roughly linearly 
with z. We therefore determine the slope of this dependence. 
In Table[TJwe compare these three terms at kiz = — 1, eval- 
uated in units of iitokiB^, as well as the value of K[/r]to. 
In Fig. |2] we show the z dependence of these three terms for 
Run B5, where R m = 68. The values of Kf as a function of 
i? m is given in the last column of TableQ] The z dependence 
of and /i f is shown in the last panel of Fig. |2] Note that 
the two profiles agree quite well. 

We point out that, near z = 0, all simulations show ei- 
ther a local reduction in the gradients of the terms on the 
RHS of Eq. (fTTT i or even a local reversal of the gradient. 
This is likely to be associated with a local reduction in dy- 
namo activity near z = 0, where kinetic helicity is zero. The 
non-uniformity of the turbulent magnetic field also leads to 
transport effects (Brandenburg & Subramanian 2005a) that 
may modify the gradient. However, we shall not pursue this 
question further here. 

Looking at Table [TJ we see that 25 • B and 2-qj ■ b 
balance each other nearly perfectly, and that only a small 
residual is then balanced by the diffusive flux divergence, 
V • J-[ . For the values of R m considered here, the terms 




-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 




Fig. 2 z dependence of the terms on the RHS of Eq. (fTTb 
in the first two panels and in Eq. ([Pfl ) for Run B5. 



2£ ■ B and 2rjj ■ b scale with i? m , while the dependence 

of V • Ti on R m is comparatively weak. If catastrophic 
quenching is to be alleviated by the magnetic helicity flux, 
one would expect that at large values of R m the terms 2£ B 

and V • J- 1 should balance. At the moment our values of 
R m are still too small by about a factor of 30-60 (assum- 
ing that the same scaling with R m persists). This result is 
compatible with that of earlier mean field models (Bran- 
denburg et al. 2009). Consequently, we see that the energy 
of the mean magnetic field decreases with increasing from 
R m = 33 to 68; see Fig. [5] For larger values of R m the 
situation is still unclear. 

In Table [TJ we also give the approximate values of 
Kf/rjto- Note that this ratio is always around 0.3 and inde- 
pendent of R m . This is the first time that an estimate for the 
diffusion coefficient of the diffusive flux has been obtained. 
There exists no theoretical prediction for value of Kf other 
than the naive expectation that such a term should be ex- 



ec) 0000 WILEY- VCH Verlag GmbH & Co. KGaA, Weinheim 



www.an-journal.org 



Astron. Nachr. / AN (0000) 



5 



10.0 



1.0 



0.1 



1.0 



1 




R' 1 




r B z /Bt q y 






. — -a 








1 1 ■ ' 1 L. 









10 



100 



0.1 



1 1 1 ■ 




: uv,/ k,u rms 


_ ~ts * - - 

Is""" 


- - A 


■ jb/k^U / 

A 






-. k,/B z eq J 

1 1 ■ ■ ■ ' 1 l_ 







10 



100 



Fig. 3 R m dependence of the normalized magnetic en- 
ergy of the mean field, (B ) / B 2 q , and the fluctuating field, 
(b 2 )/B 2 q , in the upper panel together with the normalized 
helicities of the small-scale magnetic field, a ■ bk[/B 2 q , the 
small-scale current density, j ■ b/k{B 2 q , and the small-scale 
velocity, u> ■ u/k{U 2 ms , at kiz — — 1 (i.e. in the south) in the 
lower panel. (All three helicities are negative in the north 
and positive in the south.) The shaded ares indicates that the 
solutions are different in nature, and that the simulations 
may not have run for long enough. 



pected and that its value should be of the order of r/to- This 
now allows us to state more precisely the point where the 
turbulent diffusive helicity flux becomes comparable with 
the resistive term, i.e. we assume K{ V 2 a ■ b to become com- 
parable with 2r/j ■ b. Using the relation j ■ b « kfa-b 
(Blackman & Brandenburg 2002), which is confirmed by 
the current simulations within a factor of about 2 (see the 
second panel of Fig. [3]), we find that 

K f /2 V > (kf/h) 2 , (16) 



where we have assumed that the Laplacian of a ■ b can be 
replaced by a k\ factor. Using our empirical finding, K{ w 
r?to/3, together with the definition rjto/ri « linns 
i? m /3, we arrive at the condition 

R m > 18(fc f /fci) 2 w 4600 (for m to be important), (17) 

where we have inserted the value kt jk\ = 16 for the present 
simulations. Similarly, large values of R m for alleviating 
catastrophic quenching by turbulent diffusive helicity fluxes 
were also found using mean-field modelling (Brandenburg 
et al. 2009). Unfortunately, the computing resources are still 
not sufficient to verify this in the immediate future. 



5 Gauge-dependence of helicity flux 

Let us now consider the question of gauge-dependence 

of the helicity flux. Equation ( fTTT i is obviously gauge 

— IV 

dependent. However, if, in the statistically steady state, h { 
becomes independent of time, we can average this equation 
and obtain 



■M 



= -2S ■ B - 2r]j ■ b, 



(18) 



dz 

where J^f refers to the z component of . On the RHS of 
this equation the two terms are gauge-independent. There- 
fore V • J- { must also be gauge-independent. The same 
applies also to JF^ and T^; see Eq. (0. We have con- 

— M 

firmed that, in the steady state, h f is statistically steady 
and does not show a long-term trend; cf. Fig.|4]for the three 

— M 

gauges. We note that the fluctuations of h { are typically 
much larger for the Weyl gauge than for the other two. 

We now verify the expected gauge-independence ex- 
plicitly for three different gauges: the Weyl gauge, defined 
by 

i' = 0, (19) 
the Lorenz gauge (or pseudo-Lorenz gauge defined by 

ftV = -4 V ' A ' ( 20 ) 

and the resistive gauge, defined by (TT~5T > above. We calcu- 
late the normalized magnetic helicity for both the mean and 
fluctuating parts and the respective fluxes for all the three 
gauges. These simulations are done for the Model A with 
low R m (R m *s 1.9). 

We find the transport coefficient K{ in the way described 
in the previous section. A snapshot of the mean flux T { is 
plotted in the top panel of Fig. [5] The flux is different in all 
the three gauges. However, when averaged over the horizon- 
tal directions as well as time the fluxes in the three different 
gauges agree with one another as shown in the bottom panel 
of Fig. [5] We find the transport coefficient K{ as described 
in the previous section and obtain the same value in all the 
three gauges. 



6 Conclusion 

In this paper we use a setup in which the two parts of the 
domain have different signs of kinetic and magnetic helici- 
ties. Using DNS we show that the flux of magnetic helicity 
due to small-scale fields can be described by Fickian dif- 
fusion down the gradient of this quantity. The correspond- 
ing diffusion coefficient is approximately independent of 
R m . However, in the range of R m values we have consid- 
ered here, the flux is not big enough to alleviate the catas- 
trophic quenching. The critical value of i? m for the flux to 



2 In fact, this is not the true Lorenz gauge because we use velocity of 
sound (Brandenburg & Kapyla 2007) instead of the velocity of light which 
appears in the original Lorenz gauge 
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Fig. 4 Plot of h { as a function of time in the statistically 
stationary state for k\z = — 1 (south, top panel) and k\Z = 
1 (north, bottom panel) for the three different gauges, Weyl 
gauge (open circle), Lorenz gauge (line) and resistive gauge 
(broken line). 

become important is proportional to the square of the scale 
separation ratio. In the present case, where this ratio is 16, 
the critical value of R m is estimated to be 4600. We have 
also calculated the flux and the diffusion coefficient in the 
three gauges discussed above and have found the fluxes to 
be independent of the choice of these gauges. This is ex- 
plained by the fact that in the steady state the divergence of 
magnetic helicity flux is balanced by terms that are gauge- 
independent. 

Several immediate improvements on this study spring 
to mind. One is to compare our results with the gauge- 
independent magnetic helicity of Berger & Field (1984) and 
the corresponding magnetic helicity flux. The second is to 
extend the present study to higher values of R m to under- 
stand the asymptotic behavior of the flux. Finally, it may be 
useful to compare the results for different profiles of kinetic 
helicity to see whether or not our results depend on such 
details. 
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Fig. 5 Comparison of the flux JF f (z,t) at a randomly 
chosen instant (upper panel) and its time average T f (z) for 
the three different gauges. Lorenz gauge (o), Weyl gauge (o) 
and the resistive gauge (•). The instantaneous flux is plotted 
in the top panel and the time-averaged flux is plotted in the 
bottom panel. 
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